home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Pascal / Applications / NIH Image 1.59 / 1.59 Source / Background.p < prev    next >
Encoding:
Text File  |  1995-07-24  |  39.6 KB  |  839 lines  |  [TEXT/PJMM]
  1. unit Background;
  2. {************************************************************************}
  3. {*     by Michael Castle and Janice Keller                                                                                                    *}
  4. {*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
  5. {*     e-mail address: mike.castle@umich.edu                                                                                  *}
  6. {************************************************************************}
  7. {*     Rolling ball and rolling disk algorithms inspired by Stanley Sternberg's article, "Biomedical        *}
  8. {*  Image Processing", IEEE Computer, January 1983.  Discussions with Rork Kuick and Tom               *}
  9. {*  Ford-Holevinski also contributed to the development of the algorithms and improvements in their   *}
  10. {*  efficiency.                                                                                                                                                *}
  11. {************************************************************************}
  12.  
  13. interface
  14.  
  15.     uses
  16.         Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils,Resources, Errors, Palettes, globals, Utilities, Graphics, Camera, Filters;
  17.  
  18.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  19.  
  20.  
  21. implementation
  22.  
  23.     type
  24.         IntRow = array[0..MaxLine] of Integer;
  25.         BackSubKindType = (RollingHorizontalArc, RollingVerticalArc, RollingBothArcs, RollingBall);
  26.  
  27.     var
  28.         ArcTrimPer: integer;                                  {trim off percentage of each side of the rolling ball patch}
  29.         shrinkfactor: integer;                                 {shrink the image and ball by this factor before rolling ball}
  30.         BackSubKind: BackSubKindType;                {which kind of background subtraction are we doing}
  31.         IntPlotWidth: Boolean;
  32.         Intplotwidthval: integer;
  33.         BoundRect: rect;
  34.         xxcenter, yycenter: integer;                                        {center of rectangular mask used in MinIn2DMask}
  35.         xmaskmin, ymaskmin: integer;                                   {upper left corner of mask used in AvgIn2DMask}
  36.         backgroundptr, ballptr, smallimageptr: ptr;              {ptrs to background, rolling ball, shrunk image memory}
  37.         backgroundaddr, balladdr, smallimageaddr: longint;   {addrs of background, rolling ball shrunk image memory}
  38.         patchwidth: LongInt;                                                     {x or y dimension of the rolling ball patch}
  39.         leftroll, rightroll, toproll, bottomroll: integer;         {bounds of the shrunk image}
  40.         Aborting: boolean;
  41.  
  42.  
  43.  
  44.  
  45.  
  46.     procedure RollBall;
  47. {*******************************************************************************}
  48. {*     RollBall 'rolls' a filtering object over a (shrunken) image in order to find the image's smooth continuous    *}
  49. {*  background.  For the purpose of explaining this algorithm, imagine that the 2D grayscale image has a third     *}
  50. {*  (height) dimension defined by the intensity value (0-255) at every point in the image.  The center of the      *}
  51. {*  filtering object, a patch from the top of a sphere having radius BallRadius, is moved along each scan line of     *}
  52. {*  the image so that the patch is tangent to the image at one or more points with every other point on the patch    *}
  53. {*  below the corresponding (x,y) point of the image.  Any point either on or below the patch during this process*}
  54. {*  is considered part of the background.  Shrinking the image before running this procedure is advised due to      *}
  55. {*  the fourth-degree complexity of the algorithm.  Care has been taken to avoid unnecessary operations (exp.      *}
  56. {*  multiplication inside loops) in this code.                                                                                                                *}
  57. {*******************************************************************************}
  58.         var
  59.             halfpatchwidth,                                {distance in x or y from patch center to any edge}
  60.             ptsbelowlastpatch,                           {number of points we may ignore because they were below last patch}
  61.             left, right, top, bottom,                   {}
  62.             xpt, ypt,                                           {current (x,y) point in the shrunken image}
  63.             xpt2, ypt2,                                      {current (x,y) point in the patch relative to upper left corner}
  64.             xval, yval,                                        {location in ball in shrunken image coordinates}
  65.             zdif,                                                  {difference in z (height) between point on ball and point on image}
  66.             zmin,                                                {smallest zdif for ball patch with center at current point}
  67.             zctr,                                                 {current height of the center of the sphere of which the patch is a part}
  68.             zadd: integer;                                    {height of a point on patch relative to the xy-plane of the shrunken image}
  69.             ballpt,                                               {index to chunk of memory storing the precomputed ball patch}
  70.             imgpt,                                               {index to chunk of memory storing the shrunken image}
  71.             backgrpt,                                          {index to chunk of memory storing the calculated background}
  72.             ybackgrpt,                                        {displacement to current background scan line}
  73.             ybackgrinc,                                      {distance in memory between two shrunken y-points in background}
  74.             smallimagewidth: longint;               {length of a scan line in shrunken image}
  75.             p1, p2: ptr;                                      {temporary pointers to background, ball, or small image}
  76.     begin
  77.         UpdateMeter(20, 'Finding Background...');
  78.         left := 1;
  79.         right := rightroll - leftroll - 1;
  80.         top := 1;
  81.         bottom := bottomroll - toproll - 1;
  82.         smallimagewidth := right - left + 3;
  83.         halfpatchwidth := patchwidth div 2;
  84.         ybackgrinc := shrinkfactor * (BoundRect.right - BoundRect.left);  {real dist btwn 2 adjacent (dy=1) shrunk pts}
  85.         zctr := 0;                                            {start z-center in the xy-plane}
  86.         for ypt := top to (bottom + patchwidth) do begin
  87.                 for xpt := left to (right + patchwidth) do {while patch is tangent to edges or within image...}
  88.                     begin                                           {xpt is far right edge of ball patch}
  89. {do we have to move the patch up or down to make it tangent to but not above image?...}
  90.                         zmin := 255;                              {highest could ever be 255}
  91.                         ballpt := balladdr;
  92.                         ypt2 := ypt - patchwidth;          {ypt2 is top edge of ball patch}
  93.                         imgpt := smallimageaddr + ypt2 * smallimagewidth + xpt - patchwidth;
  94.                         while ypt2 <= ypt do begin
  95.                                 xpt2 := xpt - patchwidth;      {xpt2 is far left edge of ball patch}
  96.                                 while xpt2 <= xpt do            {check every point on ball patch}
  97.                                     begin                                   {only examine points on }
  98.                                         if ((xpt2 >= left) and (xpt2 <= right) and (ypt2 >= top) and (ypt2 <= bottom)) then begin
  99.                                                 p1 := ptr(ballpt);
  100.                                                 p2 := ptr(imgpt);
  101.                                                 zdif := BAND(p2^, 255) - (zctr + BAND(p1^, 255));  {curve - circle points}
  102.                                                 if (zdif < zmin) then begin {keep most negative, since ball should always be below curve}
  103.                                                         zmin := zdif;
  104.                                                     end;
  105.                                             end;  {if xpt2,ypt2}
  106.                                         ballpt := ballpt + 1;          {step thru the ball patch memory}
  107.                                         xpt2 := xpt2 + 1;
  108.                                         imgpt := imgpt + 1;
  109.                                     end;  {while xpt2 }
  110.                                 ypt2 := ypt2 + 1;
  111.                                 imgpt := imgpt - patchwidth - 1 + smallimagewidth;
  112.                             end;  {while ypt2}
  113.                         if (zmin <> 0) then
  114.                             zctr := zctr + zmin;                {move ball up or down if we find a new minimum}
  115.                         if (zmin < 0) then
  116.                             ptsbelowlastpatch := halfpatchwidth    {ignore left half of ball patch when dz < 0}
  117.                         else
  118.                             ptsbelowlastpatch := 0;
  119. {now compare every point on ball with background,  and keep highest number}
  120.                         yval := ypt - patchwidth;
  121.                         ypt2 := 0;
  122.                         ballpt := balladdr;
  123.                         ybackgrpt := backgroundaddr + (yval - top + 1) * ybackgrinc;
  124.                         while ypt2 <= patchwidth do begin
  125.                                 xval := xpt - patchwidth + ptsbelowlastpatch;
  126.                                 xpt2 := ptsbelowlastpatch;
  127.                                 ballpt := ballpt + ptsbelowlastpatch;
  128.                                 backgrpt := ybackgrpt + (xval - left + 1) * shrinkfactor;
  129.                                 while xpt2 <= patchwidth do begin     {for all the points in the ball patch}
  130.                                         if ((xval >= left) and (xval <= right) and (yval >= top) and (yval <= bottom)) then begin
  131.                                                 p1 := ptr(ballpt);
  132.                                                 zadd := zctr + BAND(p1^, 255);
  133.                                                 p1 := ptr(backgrpt);
  134.                                                 if (zadd > BAND(p1^, 255)) then  {keep largest adjustment}
  135.                                                     p1^ := zadd;
  136.                                             end;
  137.                                         ballpt := ballpt + 1;
  138.                                         xval := xval + 1;
  139.                                         xpt2 := xpt2 + 1;
  140.                                         backgrpt := backgrpt + shrinkfactor;     {move to next point in x}
  141.                                     end;  {while xpt2}
  142.                                 yval := yval + 1;
  143.                                 ypt2 := ypt2 + 1;
  144.                                 ybackgrpt := ybackgrpt + ybackgrinc;       {move to next point in y}
  145.                             end;  {while ypt2}
  146.                     end;  {for xpt }
  147.                 if ((ypt mod 5) = 0) or not FasterBackgroundSubtraction then begin
  148.                         UpdateMeter(20 + (ord4(ypt - top) * 70) div (bottom + patchwidth - top), 'Finding Background...');
  149.                         if CommandPeriod then begin
  150.                                 beep;
  151.                                 Aborting := true;
  152.                                 Exit(RollBall);
  153.                             end;
  154.                     end;
  155.             end;  {for ypt}
  156.     end;
  157.  
  158.  
  159.     function MinIn2DMask {(xmaskmin,ymaskmin: integer)}
  160.         : integer;
  161. {*******************************************************************************}
  162. {*     MinInMask finds the minimum pixel value in a shrinkfactor X shrinkfactor mask.                                           *}
  163. {*******************************************************************************}
  164.         var
  165.             i, j,                                           {loop indices to step through mask}
  166.             thispixel,                                  {value at current pixel in mask}
  167.             min,                                          {temporary minimum value in mask}
  168.             nextrowoffset: integer;             {distance in memory from end of mask in this row to beginning in next}
  169.             paddr: longint;                           {address of current mask pixel}
  170.             p: ptr;                                        {pointer to current pixel in mask}
  171.     begin
  172.         with info^ do begin
  173.                 min := 255;
  174.                 nextrowoffset := bytesperrow - shrinkfactor;
  175.                 paddr := ord4(PicBaseAddr) + ymaskmin * bytesperrow + xmaskmin;
  176.                 for j := 1 to shrinkfactor do begin
  177.                         for i := 1 to shrinkfactor do begin
  178.                                 p := ptr(paddr);
  179.                                 thispixel := BAND(p^, 255);
  180.                                 if (thispixel < min) then
  181.                                     min := thispixel;
  182.                                 paddr := paddr + 1;
  183.                             end;     {for i}
  184.                         paddr := paddr + nextrowoffset;
  185.                     end;     {for j}
  186.                 MinIn2DMask := min;
  187.             end; {with}
  188.     end;
  189.  
  190.  
  191.     procedure GetRollingBall;
  192. {******************************************************************************}
  193. {*     This procedure computes the location of each point on the rolling ball patch relative to the center of the     *}
  194. {*  sphere containing it.  The patch is located in the top half of this sphere.  The vertical axis of the sphere         *}
  195. {*  passes through the center of the patch.  The projection of the patch in the xy-plane below is a square.           *}
  196. {******************************************************************************}
  197.         var
  198.             rsquare,                                                                         {rolling ball radius squared}
  199.             xtrim,                                                                            {# of pixels trimmed off each end of ball to make patch}
  200.             xval, yval,                                                                     {x,y-values on patch relative to center of rolling ball}
  201.             smallballradius, diam,                                                  {radius and diameter of rolling ball}
  202.             temp,                                                                             {value must be >=0 to take square root}
  203.             halfpatchwidth: integer;                                                {distance in x or y from center of patch to any edge}
  204.             i,                                                                                    {index into rolling ball patch memory}
  205.             ballsize: LongInt;                                                           {size of rolling ball memory}
  206.             p: ptr;                                                                            {pointer to current point on the ball patch}
  207.     begin
  208.         smallballradius := ballradius div shrinkfactor;           {operate on small-sized image with small-sized ball}
  209.         if smallballradius < 1 then
  210.             smallballradius := 1;
  211.         rsquare := smallballradius * smallballradius;
  212.         diam := smallballradius * 2;
  213.         xtrim := (ArcTrimPer * diam) div 100;                      {only use a patch of the rolling ball}
  214.         patchwidth := diam - xtrim - xtrim;
  215.         halfpatchwidth := smallballradius - xtrim;                   {this is half the patch width}
  216.         ballsize := (patchwidth + 1) * (patchwidth + 1);
  217.         ballptr := NewPtrClear(ballsize);
  218.         p := ballptr;
  219.         for i := 0 to ballsize - 1 do begin
  220.                 xval := i mod (patchwidth + 1) - halfpatchwidth;
  221.                 yval := i div (patchwidth + 1) - halfpatchwidth;
  222.                 temp := rsquare - (xval * xval) - (yval * yval);
  223.                 if (temp >= 0) then
  224.                     p^ := round(sqrt(temp))
  225.                 else
  226.                     p^ := 0;
  227.                 p := ptr(ord4(p) + 1);
  228.             end;
  229.     end;
  230.  
  231.  
  232.     procedure InterpolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  233. {******************************************************************************}
  234. {*     This procedure uses bilinear interpolation to find the points in the full-scale background given the points *}
  235. {*  from the shrunken image background.  Since the shrunken background is found from an image composed of    *}
  236. {*  minima (over a sufficiently large mask), it is certain that no point in the full-scale interpolated                 *}
  237. {*  background has a higher pixel value than the corresponding point in the original image.                                  *}
  238. {******************************************************************************}
  239.         var
  240.             i, ii,                                                   {horizontal loop indices}
  241.             j, jj,                                                  {vertical loop indices}
  242.             hloc, vloc,                                          {position of current pixel in calculated background}
  243.             vinc,                                                   {memory offset from current calculated pos to current interpolated pos}
  244.             lastvalue, nextvalue: integer;           {calculated pixel values between which we are interpolating}
  245.             p,                                                        {pointer to current interpolated pixel value}
  246.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values between which we are interpolating}
  247.             width: LongInt;
  248.     begin
  249.         vloc := 0;
  250.         with BoundRect do begin
  251.                 width := right - left;
  252.                 for j := 1 to bottomroll - toproll - 1 do begin     {interpolate to find background interior}
  253.                         hloc := 0;
  254.                         vloc := vloc + shrinkfactor;
  255.                         for i := 1 to rightroll - leftroll do begin
  256.                                 hloc := hloc + shrinkfactor;
  257.                                 bgnextptr := ptr(backgroundaddr + vloc * width + hloc);
  258.                                 bglastptr := ptr(ord4(bgnextptr) - shrinkfactor);
  259.                                 nextvalue := BAND(bgnextptr^, 255);
  260.                                 lastvalue := BAND(bglastptr^, 255);
  261.                                 for ii := 1 to shrinkfactor - 1 do begin     {interpolate horizontally}
  262.                                         p := ptr(ord4(bgnextptr) - ii);
  263.                                         p^ := lastvalue + (shrinkfactor - ii) * (nextvalue - lastvalue) div shrinkfactor;
  264.                                     end;     {for ii}
  265.                                 for ii := 0 to shrinkfactor - 1 do begin     {interpolate vertically}
  266.                                         bglastptr := ptr(backgroundaddr + (vloc - shrinkfactor) * width + hloc - ii);
  267.                                         bgnextptr := ptr(backgroundaddr + vloc * width + hloc - ii);
  268.                                         lastvalue := BAND(bglastptr^, 255);
  269.                                         nextvalue := BAND(bgnextptr^, 255);
  270.                                         vinc := 0;
  271.                                         for jj := 1 to shrinkfactor - 1 do begin
  272.                                                 vinc := vinc - right + left;
  273.                                                 p := ptr(ord4(bgnextptr) + vinc);
  274.                                                 p^ := lastvalue + (shrinkfactor - jj) * (nextvalue - lastvalue) div shrinkfactor;
  275.                                             end;     {for jj}
  276.                                     end;     {for ii}
  277.                             end;     {for i}
  278.                     end;     {for j}
  279.             end;   {with boundrect}
  280.     end;
  281.  
  282.  
  283.     procedure ExtrapolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  284. {******************************************************************************}
  285. {*     This procedure uses linear extrapolation to find pixel values on the top, left, right, and bottom edges of      *}
  286. {*  the background.  First it finds the top and bottom edge points by extrapolating from the edges of the                *}
  287. {*  calculated and interpolated background interior.  Then it uses the edge points on the new calculated,               *}
  288. {*  interpolated, and extrapolated background to find all of the left and right edge points.  If extrapolation yields *}
  289. {*  values below zero or above 255, then they are set to zero and 255 respectively.                                              *}
  290. {******************************************************************************}
  291.         var
  292.             ii, jj,                                                 {horizontal and vertical loop indices}
  293.             hloc, vloc,                                          {position of current pixel in calculated/interpolated background}
  294.             edgeslope,                                          {difference of last two consecutive pixel values on an edge}
  295.             pvalue,                                               {current extrapolated pixel value}
  296.             lastvalue, nextvalue: integer;           {calculated pixel values from which we are extrapolating}
  297.             p,                                                        {pointer to current extrapolated pixel value}
  298.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values from which we are extrapolating}
  299.             width: LongInt;
  300.     begin
  301.         with BoundRect do begin
  302.                 width := right - left;
  303.                 for hloc := shrinkfactor to shrinkfactor * (rightroll - leftroll) - 1 do begin     {extrapolate on top and bottom}
  304.                         bglastptr := ptr(backgroundaddr + shrinkfactor * width + hloc);
  305.                         bgnextptr := ptr(backgroundaddr + (shrinkfactor + 1) * width + hloc);
  306.                         lastvalue := BAND(bglastptr^, 255);
  307.                         nextvalue := BAND(bgnextptr^, 255);
  308.                         edgeslope := nextvalue - lastvalue;
  309.                         p := bglastptr;
  310.                         pvalue := lastvalue;
  311.                         for jj := 1 to shrinkfactor do begin
  312.                                 p := ptr(ord4(p) - right + left);
  313.                                 pvalue := pvalue - edgeslope;
  314.                                 if (pvalue < 0) then
  315.                                     p^ := 0
  316.                                 else if (pvalue > 255) then
  317.                                     p^ := 255
  318.                                 else
  319.                                     p^ := pvalue;
  320.                             end;     {for jj}
  321.                         bglastptr := ptr(backgroundaddr + (shrinkfactor * (bottomroll - toproll - 1) - 1) * width + hloc);
  322.                         bgnextptr := ptr(backgroundaddr + shrinkfactor * (bottomroll - toproll - 1) * width + hloc);
  323.                         lastvalue := BAND(bglastptr^, 255);
  324.                         nextvalue := BAND(bgnextptr^, 255);
  325.                         edgeslope := nextvalue - lastvalue;
  326.                         p := bgnextptr;
  327.                         pvalue := nextvalue;
  328.                         for jj := 1 to (bottom - top - 1) - shrinkfactor * (bottomroll - toproll - 1) do begin
  329.                                 p := ptr(ord4(p) + right - left);
  330.                                 pvalue := pvalue + edgeslope;
  331.                                 if (pvalue < 0) then
  332.                                     p^ := 0
  333.                                 else if (pvalue > 255) then
  334.                                     p^ := 255
  335.                                 else
  336.                                     p^ := pvalue;
  337.                             end;     {for jj}
  338.                     end;     {for hloc}
  339.                 for vloc := top to bottom - 1 do begin     {extrapolate on left and right}
  340.                         bglastptr := ptr(backgroundaddr + (vloc - top) * width + shrinkfactor);
  341.                         bgnextptr := ptr(ord4(bglastptr) + 1);
  342.                         lastvalue := BAND(bglastptr^, 255);
  343.                         nextvalue := BAND(bgnextptr^, 255);
  344.                         edgeslope := nextvalue - lastvalue;
  345.                         p := bglastptr;
  346.                         pvalue := lastvalue;
  347.                         for ii := 1 to shrinkfactor do begin
  348.                                 p := ptr(ord4(p) - 1);
  349.                                 pvalue := pvalue - edgeslope;
  350.                                 if (pvalue < 0) then
  351.                                     p^ := 0
  352.                                 else if (pvalue > 255) then
  353.                                     p^ := 255
  354.                                 else
  355.                                     p^ := pvalue;
  356.                             end;     {for ii}
  357.                         bgnextptr := ptr(backgroundaddr + (vloc - top) * width + shrinkfactor * (rightroll - leftroll - 1) - 1);
  358.                         bglastptr := ptr(ord4(bgnextptr) - 1);
  359.                         lastvalue := BAND(bglastptr^, 255);
  360.                         nextvalue := BAND(bgnextptr^, 255);
  361.                         edgeslope := nextvalue - lastvalue;
  362.                         p := bgnextptr;
  363.                         pvalue := nextvalue;
  364.                         for ii := 1 to (right - left - 1) - shrinkfactor * (rightroll - leftroll - 1) + 1 do begin
  365.                                 p := ptr(ord4(p) + 1);
  366.                                 pvalue := pvalue + edgeslope;
  367.                                 if (pvalue < 0) then
  368.                                     p^ := 0
  369.                                 else if (pvalue > 255) then
  370.                                     p^ := 255
  371.                                 else
  372.                                     p^ := pvalue;
  373.                             end;     {for ii}
  374.                     end;     {for vloc}
  375.             end;   {with BoundRect}
  376.     end;
  377.  
  378.  
  379.     procedure SubtractBackground2D;
  380. {*****************************************************************************}
  381. {*     This procedure subtracts each pixel from the calculated/interpolated/extrapolated background from the  *}
  382. {*  corresponding pixel value in the original image.  The resulting image is stored in place of the original        *}
  383. {*  image.  Any pixel subtractions with results below zero are given the value zero.                                           *}
  384. {*****************************************************************************}
  385.         var
  386.             hloc, vloc,                                          {current pixel location in image and background}
  387.             pvalue: integer;                                 {difference at current pixel location}
  388.             offset,                                                 {offset in memory from beginning of original image to current scan line}
  389.             backgrpt: LongInt;                              {offset to current point in background}
  390.             p: ptr;                                                {temporary pointer to image or background points}
  391.             Databand: Linetype;                           {current scan line in image}
  392.             ControlKey: boolean;
  393.     begin
  394.         backgrpt := 0;
  395.         ControlKey := ControlKeyDown;
  396.         with Info^, BoundRect do begin
  397.                 for vloc := top to bottom - 1 do begin
  398.                         GetLine(0, vloc, pixelsperline, Databand);
  399.                         for hloc := left to right - 1 do begin
  400.                                 p := ptr(backgroundaddr + backgrpt);
  401.                                 pvalue := Databand[hloc] - BAND(p^, 255);
  402.                                 if ControlKey then
  403.                                     pvalue := BAND(p^, 255);
  404.                                 if pvalue < 0 then
  405.                                     Databand[hloc] := 0
  406.                                 else
  407.                                     Databand[hloc] := pvalue;
  408.                                 backgrpt := backgrpt + 1;
  409.                             end;     {for}
  410.                         offset := vloc * BytesPerRow;
  411.                         p := ptr(ord4(PicBaseAddr) + offset);
  412.                         BlockMove(@Databand, p, pixelsperline);
  413.                     end;  {for}
  414.             end;     {with}
  415.     end;
  416.  
  417.  
  418.     procedure Background2D;
  419. {******************************************************************************}
  420. {*     This procedure implements a rolling-ball algorithm for the removal of smooth continuous background       *}
  421. {*  from a two-dimensional gel image.  It rolls the ball (actually a square patch on the top of a sphere) on a       *}
  422. {*  low-resolution (by a factor of 'shrinkfactor' times) copy of the original image in order to increase speed     *}
  423. {*  with little loss in accuracy.  It uses interpolation and extrapolation to blow the shrunk image to full size.     *}
  424. {******************************************************************************}
  425.         var
  426.             tport: Grafptr;
  427.             i,                                     {loop index for shrunk image memory}
  428.             backgroundsize,              {size of the background memory}
  429.             smallimagesize: LongInt;     {size of the shrunk image memory}
  430.             p: ptr;                             {pointer to current pixel in shrunk image memory}
  431.             table: FateTable;             {not used}
  432.             width: LongInt;
  433.     begin
  434.         ShowWatch;
  435.         UpdateMeter(0, 'Building Rolling Ball...');
  436.         GetPort(tPort);
  437.         with Info^ do begin
  438.                 SetPort(GrafPtr(osPort));
  439.                 BoundRect := roiRect;
  440.             end;
  441.         GetRollingBall;                                                                  {precompute the rolling ball}
  442.         UpdateMeter(3, 'Building Rolling Ball...');
  443.         balladdr := ord4(ballptr);
  444.         with BoundRect do begin
  445.                 width := right - left;
  446.                 leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  447.                 rightroll := right div shrinkfactor - 1;                      {on which to roll ball}
  448.                 toproll := top div shrinkfactor;
  449.                 bottomroll := bottom div shrinkfactor - 1;
  450.                 backgroundsize := width * (bottom - top);
  451.                 backgroundptr := NewPtrClear(backgroundsize);
  452.                 Aborting := backgroundptr = nil;
  453.                 backgroundaddr := ord4(backgroundptr);
  454.                 smallimagesize := ord4(rightroll - leftroll + 1);
  455.                 smallimagesize := smallimagesize * (bottomroll - toproll + 1);
  456.                 smallimageptr := NewPtrClear(smallimagesize);
  457.                 Aborting := Aborting or (smallimageptr = nil);
  458.                 smallimageaddr := ord4(smallimageptr);
  459.                 if not aborting then begin
  460.                         UpdateMeter(6, 'Smoothing Image ');
  461.                         filter(unweightedAvg, 1, table);                                {smooth image before shrinking}
  462.                         UpdateMeter(10, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  463.                         for i := 0 to smallimagesize - 1 do begin                {create a lower resolution image for ball-rolling}
  464.                                 p := ptr(smallimageaddr + i);
  465.                                 xmaskmin := left + shrinkfactor * (i mod (rightroll - leftroll + 1));
  466.                                 ymaskmin := top + shrinkfactor * (i div (rightroll - leftroll + 1));
  467.                                 p^ := MinIn2DMask;                            {each point in small image is minimum of its neighborhood}
  468.                             end;
  469.                         if not aborting then begin
  470.                                 Undo;        {restore original unsmoothed image}
  471.                                 RollBall;
  472.                             end;
  473.                     end
  474.                 else
  475.                     beep;
  476.                 if not Aborting then begin
  477.                         UpdateMeter(90, 'Interpolating Background...');
  478.                         InterpolateBackground2D;                              {interpolate to find background interior}
  479.                         UpdateMeter(95, 'Extrapolating Background...');
  480.                         ExtrapolateBackground2D;                             {extrapolate on top and bottom}
  481.                         UpdateMeter(98, 'Subtracting Background...');
  482.                         SubtractBackground2D;                                  {subtract background from original image}
  483.                         UpdateMeter(100, 'Subtracting Background...');
  484.                     end;
  485.             end;   {with boundrect}
  486.         DisposePtr(backgroundptr);                           {free up background, rolling ball, shrunk image memory}
  487.         DisposePtr(ballptr);
  488.         DisposePtr(smallimageptr);
  489.         DisposeWindow(MeterWindow);
  490.         MeterWindow := nil;
  491.         SetPort(tPort);
  492.     end;
  493.  
  494.  
  495.     procedure RollArc (left, rightminusone, diam: integer; var background, ballpoints: IntRow; var Dataline: Linetype);
  496.         var
  497.             xpt, xpt2, xval, ydif, ymin, yctr, bpt, yadd: integer;
  498.     begin
  499.         for xpt := left to rightminusone do begin
  500.                 background[xpt] := -255;         {init background curve to minimum values}
  501.             end;
  502.         yctr := 0;                                   {start y-center at the x axis}
  503.         for xpt := left to (rightminusone + diam - 1) do {while semicircle is tangent to edges or within curve...}
  504.             begin                                       {xpt is far right edge of semi-circle}
  505. {do we have to move the circle?...}
  506.                 ymin := 256;                          {highest could ever be 255}
  507.                 bpt := 0;
  508.                 xpt2 := xpt - diam;                {xpt2 is far left edge of semi-circle}
  509.                 while xpt2 <= xpt do            {check every point on semicircle}
  510.                     begin
  511.                         if ((xpt2 >= left) and (xpt2 <= rightminusone)) then begin  {only examine points on curve}
  512.                                 ydif := dataline[xpt2] - (yctr + ballpoints[bpt]);                {curve minus circle points}
  513.                                 if (ydif < ymin) then begin {keep most negative, since ball should always be below curve}
  514.                                         ymin := ydif;
  515.                                     end;
  516.                             end;  {if xpt2 }
  517.                         bpt := bpt + 1;
  518.                         xpt2 := xpt2 + 1;
  519.                     end;  {while xpt2 }
  520.                 if (ymin <> 256) then{if we found a new minimum...}
  521.                     yctr := yctr + ymin;   {move circle up or down}
  522. {now compare every point on semi with background,  and keep highest number}
  523.                 xval := xpt - diam;
  524.                 xpt2 := 0;
  525.                 while xpt2 <= diam do begin  {for all the points in the semicircle}
  526.                         if ((xval >= left) and (xval <= rightminusone)) then begin
  527.                                 yadd := yctr + ballpoints[xpt2];
  528.                                 if (yadd > background[xval]) then  {keep largest adjustment}
  529.                                     background[xval] := yadd;
  530.                             end;
  531.                         xval := xval + 1;
  532.                         xpt2 := xpt2 + 1;
  533.                     end;  {while xpt2}
  534.             end;  {for xpt }
  535.     end;
  536.  
  537.  
  538.     function MinIn1DMask (var Databand: LineType; xcenter: integer): integer;
  539. {*******************************************************************************}
  540. {*     MinIn1DMask finds the minimum pixel value in a (2*shrinkfactor-1) mask about the point xcenter in the *}
  541. {*  current line.  This code must run FAST because it gets called OFTEN!                                                                   *}
  542. {*******************************************************************************}
  543.         var
  544.             i,                                                                              {index to pixels in the mask}
  545.             temp: integer;                                                          {temporary minimum value}
  546.     begin
  547.         temp := Databand[xcenter - shrinkfactor + 1];
  548.         for i := xcenter - shrinkfactor + 2 to xcenter + shrinkfactor - 1 do
  549.             if (Databand[i] < temp) then
  550.                 temp := Databand[i];
  551.         MinIn1DMask := temp;
  552.     end;
  553.  
  554.  
  555. {******************************************************************************}
  556. {*     This procedure computes the location of each point on the rolling arc relative to the center of the circle     *}
  557. {*  containing it.  The arc is located in the top half of this circle.  The vertical axis of the circle passes through  *}
  558. {*  the midpoint of the arc.  The projection of the arc on the x-axis below is a line segment.                                 *}
  559. {******************************************************************************}
  560.     procedure GetRollingArc (var arcpoints: IntRow; var arcwidth: integer);
  561.         var
  562.             xpt,                                                                                 {x-point along arc}
  563.             xval,                                                                               {x-point in arc array}
  564.             rsquare,                                                                         {shrunken arc radius squared}
  565.             xtrim,                                                                            {points to be trimmed off each end of arc}
  566.             smallballradius,                                                            {radius of shrunken arc which actually rolls}
  567.             diam: integer;                                                                 {diameter of shrunken arc's circle}
  568.     begin
  569.         smallballradius := ballradius div shrinkfactor;            { operate on small-sized image with small-sized ball}
  570.         rsquare := smallballradius * smallballradius;
  571.         for xpt := -smallballradius to smallballradius do        { find the ballpoints for arc based at  (x,y)=(0,0) }
  572.             begin
  573.                 xval := xpt + smallballradius;                                     {offset, can't have negative index}
  574.                 arcpoints[xval] := round(sqrt(abs(rsquare - (xpt * xpt))));  {Ys are positive, top half of circle}
  575.             end;
  576.         diam := smallballradius * 2;
  577.         xtrim := (ArcTrimPer * diam) div 100;                       {how many points to trim off each end}
  578.         arcwidth := diam - xtrim - xtrim;
  579.         for xpt := -smallballradius to smallballradius - xtrim - xtrim do begin
  580.                 xval := xpt + smallballradius;
  581.                 arcpoints[xval] := arcpoints[xval + xtrim];
  582.             end;
  583.         for xpt := smallballradius - xtrim - xtrim + 1 to smallballradius do begin
  584.                 xval := xpt + smallballradius;
  585.                 arcpoints[xval] := 0;
  586.             end;
  587.     end;
  588.  
  589.  
  590.     procedure ExtrapolateBackground1D (var Backline, Dataline: LineType; background: IntRow; leftroll, rightroll: integer);
  591. {******************************************************************************}
  592. {*     This procedure uses linear extrapolation to find pixel values on the left and right edges of the current        *}
  593. {*  line of the background.  It finds the edge points by extrapolating from the edges of the calculated and               *}
  594. {*  interpolated background interior.  If extrapolation yields values below zero or above 255, then they are set *}
  595. {*  to zero and 255 respectively.                                                                                                                               *}
  596. {******************************************************************************}
  597.         var
  598.             i,                                                                             {index to edges of background array}
  599.             hloc,                                                                       {}
  600.             edgeslope: integer;                                                 {}
  601.     begin
  602.         with BoundRect do begin
  603.                 edgeslope := (background[leftroll + 1] - background[leftroll + 2]) div shrinkfactor;
  604.                 for i := left to shrinkfactor * (leftroll + 1) - 1 do begin     {extrapolate on left edge}
  605.                         hloc := shrinkfactor * (leftroll + 1) - 1 + left - i;
  606.                         if (Backline[hloc + 1] + edgeslope < 0) then
  607.                             Backline[hloc] := 0
  608.                         else if (Backline[hloc + 1] + edgeslope > Dataline[hloc]) then
  609.                             Backline[hloc] := Dataline[hloc]
  610.                         else
  611.                             Backline[hloc] := Backline[hloc + 1] + edgeslope;
  612.                     end;
  613.                 edgeslope := (background[rightroll] - background[rightroll - 1]) div shrinkfactor;
  614.                 for hloc := shrinkfactor * rightroll to right - 1 do begin     {extrapolate on right edge}
  615.                         if (Backline[hloc - 1] + edgeslope < 0) then
  616.                             Backline[hloc] := 0
  617.                         else if (Backline[hloc - 1] + edgeslope > Dataline[hloc]) then
  618.                             Backline[hloc] := Dataline[hloc]
  619.                         else
  620.                             Backline[hloc] := Backline[hloc - 1] + edgeslope;
  621.                     end;
  622.             end;     {with}
  623.     end;
  624.  
  625.     procedure Background1D;
  626.         var
  627.             tport: Grafptr;
  628.             hloc, arcwidth, leftroll, rightroll, numpixels: integer;
  629.             left, right, top, bottom: integer;                      {image bounds; ROTATED if RollingVerticalArc}
  630.             i, j, maskwidth: integer;
  631.             background, arcpoints: IntRow;
  632.             vloc, offset: LongInt;
  633.             p: ptr;
  634.             Dataline, Backline, Smalldataline: Linetype;
  635.             str: str255;
  636.     begin
  637.         ShowWatch;
  638.         UpdateMeter(0, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  639.         GetPort(tPort);
  640.         with Info^ do begin
  641.                 SetPort(GrafPtr(osPort));
  642.                 BoundRect := roiRect;
  643.             end;
  644.         GetRollingArc(arcpoints, arcwidth);
  645.         maskwidth := shrinkfactor + shrinkfactor - 1;
  646.         case BackSubKind of
  647.             RollingHorizontalArc:  begin
  648.                     left := BoundRect.left;
  649.                     top := BoundRect.top;
  650.                     right := BoundRect.right;
  651.                     bottom := BoundRect.bottom;
  652.                     numpixels := Info^.pixelsperline;
  653.                     str := 'Rolling Disk Horizontally...';
  654.                 end;
  655.             RollingVerticalArc:  begin
  656.                     left := BoundRect.top;
  657.                     top := BoundRect.left;
  658.                     right := BoundRect.bottom;
  659.                     bottom := BoundRect.right;
  660.                     numpixels := Info^.nlines;
  661.                     str := 'Rolling Disk Vertically...';
  662.                 end;
  663.         end;     {case}
  664.         leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  665.         rightroll := right div shrinkfactor - 1;                      {on which to roll arc}
  666.         for vloc := top to bottom - 1 do begin  {for ROI}
  667.                 case BackSubKind of
  668.                     RollingHorizontalArc: 
  669.                         GetLine(0, vloc, numpixels, Dataline);
  670.                     RollingVerticalArc: 
  671.                         GetColumn(vloc, 0, numpixels, Dataline);
  672.                 end;     {case}
  673.                 for i := leftroll + 1 to rightroll do
  674.                     smalldataline[i] := MinIn1DMask(Dataline, shrinkfactor * i - 1);
  675.                 RollArc(leftroll + 1, rightroll, arcwidth, background, arcpoints, smalldataline);  {roll arc on one line}
  676.                 for i := leftroll + 1 to rightroll do begin           {interpolate to find interior background points}
  677.                         hloc := shrinkfactor * i - 1;
  678.                         Backline[hloc] := background[i];
  679.                         for j := 1 to shrinkfactor - 1 do
  680.                             Backline[hloc - j] := background[i - 1] + (shrinkfactor - j) * (background[i] - background[i - 1]) div shrinkfactor;
  681.                     end;
  682.                 ExtrapolateBackground1D(Backline, Dataline, background, leftroll, rightroll);
  683.                 for i := left to right - 1 do begin                                {subtract background from current scan line}
  684.                         Dataline[i] := Dataline[i] - Backline[i];
  685.                         if Dataline[i] < 0 then
  686.                             Dataline[i] := 0;
  687.                     end;
  688.                 case BackSubKind of
  689.                     RollingHorizontalArc: 
  690.                         with Info^ do begin
  691.                                 offset := vloc * BytesPerRow;
  692.                                 p := ptr(ord4(PicBaseAddr) + offset);
  693.                                 BlockMove(@Dataline, p, numpixels);            {fast whole line write}
  694.                             end;
  695.                     RollingVerticalArc: 
  696.                         PutColumn(vloc, 0, numpixels, Dataline);         {slow whole column write}
  697.                 end;     {case}
  698.                 if ((vloc mod 8) = 0) and (vloc > 16) then begin
  699.                         UpdateMeter(((vloc - top) * 100) div (bottom - top - 1), str);
  700.                         if CommandPeriod then begin
  701.                                 beep;
  702.                                 Aborting := true;
  703.                                 leave;
  704.                             end;
  705.                     end;
  706.             end;
  707.         UpdateMeter(100, str);
  708.         DisposeWindow(MeterWindow);
  709.         MeterWindow := nil;
  710.         SetPort(tPort);
  711.     end;
  712.  
  713.     procedure SetUpGel;
  714.         var
  715.             frame: rect;
  716.             AutoSelectAll: boolean;
  717.             p: ptr;
  718.     begin
  719.         if NotinBounds or NotRectangular then
  720.             exit(SetUpGel);
  721.         StopDigitizing;
  722.         AutoSelectAll := not Info^.RoiShowing;
  723.         if AutoSelectAll then
  724.             SelectAll(false);
  725.         SetupUndoFromClip;
  726.         with info^ do begin
  727.                 frame := roiRect;
  728.                 if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
  729.                     ApplyLookupTable;
  730.                 changes := true;
  731.             end;
  732.         case BackSubKind of
  733.             RollingHorizontalArc, RollingVerticalArc: 
  734.                 Background1D;    {--------------> call background subtract <-------------------}
  735.             RollingBall: 
  736.                 Background2D;
  737.             RollingBothArcs:  begin
  738.                     BackSubKind := RollingHorizontalArc;           {remove horizontal streaks}
  739.                     Background1D;
  740.                     BackSubKind := RollingVerticalArc;               {remove vertical streaks}
  741.                     if not aborting then
  742.                         Background1D;
  743.                     BackSubKind := RollingBothArcs;                   {leave BackSubKind as we found it}
  744.                 end;
  745.         end;     {case}
  746.         UpdatePicWindow;
  747.         SetUpRoiRect;
  748.         WhatToUndo := UndoFilter;
  749.         Info^.changes := true;
  750.         ShowWatch;
  751.         if AutoSelectAll then
  752.             KillRoi;
  753.         if Aborting then begin
  754.                 Undo;
  755.                 WhatToUndo := NothingToUndo;
  756.                 UpdatePicWindow;
  757.             end;
  758.     end;
  759.  
  760.  
  761.     procedure GetBallRadius;
  762.         var
  763.             SaveRadius: integer;
  764.             canceled: boolean;
  765.     begin
  766.         SaveRadius := BallRadius;
  767.         BallRadius := GetInt('Rolling BallRadius:', BallRadius, canceled);
  768.         if (BallRadius < 1) or (BallRadius > 319) or canceled then begin
  769.                 BallRadius := SaveRadius;
  770.                 if not canceled then
  771.                     beep;
  772.             end;
  773.     end;
  774.  
  775.  
  776.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  777.         var
  778.             map_array: Ptr;
  779.     begin
  780.         if MenuItem <= RemoveStreaksItem then
  781.             if not CheckCalibration
  782.                 then exit(DoBackgroundMenuEvent);
  783.         MeterWindow := nil;
  784.         Aborting := false;
  785.         case MenuItem of
  786.             HorizontalItem, VerticalItem:  begin
  787.                     if FasterBackgroundSubtraction then begin
  788.                             ArcTrimPer := 20;
  789.                             shrinkfactor := 4;
  790.                         end
  791.                     else begin
  792.                             ArcTrimPer := 10;
  793.                             shrinkfactor := 2;
  794.                         end;
  795.                     if MenuItem = HorizontalItem then
  796.                         BackSubKind := RollingHorizontalArc
  797.                     else
  798.                         BackSubKind := RollingVerticalArc;
  799.                     SetUpGel;
  800.                 end;
  801.             Sub2DItem:  begin
  802.                     if FasterBackgroundSubtraction then begin
  803.                             if BallRadius > 15 then begin
  804.                                     ArcTrimPer := 20;     {trim 40% in x and y}
  805.                                     shrinkfactor := 8;
  806.                                 end
  807.                             else begin
  808.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  809.                                     shrinkfactor := 4;
  810.                                 end
  811.                         end
  812.                     else begin  {faster not checked}
  813.                             if BallRadius > 15 then begin
  814.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  815.                                     shrinkfactor := 4;
  816.                                 end
  817.                             else begin
  818.                                     ArcTrimPer := 12;   {trim 24% in x and y}
  819.                                     ShrinkFactor := 2;
  820.                                 end
  821.                         end;
  822.                     BackSubKind := RollingBall;
  823.                     SetUpGel;
  824.                 end;
  825.             RemoveStreaksItem:  begin
  826.                     ArcTrimPer := 20;
  827.                     shrinkfactor := 4;
  828.                     BackSubKind := RollingBothArcs;
  829.                     SetUpGel;
  830.                 end;
  831.             FasterItem: 
  832.                 FasterBackgroundSubtraction := not FasterBackgroundSubtraction;
  833.             RadiusItem: 
  834.                 GetBallRadius;
  835.         end; {case}
  836.     end;
  837.  
  838.  
  839. end.